%Tx = plist(1);
%Ty = plist(2);

%C =  plist(4); %%C*v_temp = v_max
%z_min = plist(5);
%z_max = plist(6);
%dz = plist(7);
%rho = plist(8);
%E0=plist(9);
%E_zwidth=plist(10);
%E_zpos=plist(11);
%E_tau=plist(12);
%E_twidth=plist(13);
%E_tpos=plist(14);
%B0 = plist(15);
%t_start = plist(16);
%t_end = plistt(17);

function answer = stupid_kinetic(plist)
global PCOPY
global PCONSTANTS
PCOPY = plist;


e = 4.8032043e-10;
m = 9.1093822e-28;
c = 3e10;

B0 = plist(15);
E0=plist(9);
w_ce = e*B0/(m*c);
z_min = plist(5);
z_max = plist(6);
dz = plist(7);
nz = ceil(z_max-z_min)/dz;


PCONSTANTS(1)=e;
PCONSTANTS(2)=m;
PCONSTANTS(3)=c;
PCONSTANTS(4)=w_ce;
PCONSTANTS(5)=nz;


E_zwidth = plist(10);
E_tau_reduced=w_ce*plist(12);
E_twidth_reduced=w_ce*plist(13);
E_tpos_reduced=w_ce*plist(14);
PCOPY(12) = E_tau_reduced;
PCOPY(13) = E_twidth_reduced;
PCOPY(14) = E_tpos_reduced;

E_reduced = e*E0/(m*w_ce);
PCOPY(9) = E_reduced;
[f0, ~] = generate_initial_distr(PCOPY);

dt = 2*pi/w_ce;

PCONSTANTS(6)=dt;
PCONSTANTS(7) = 1/1.6e-12;
erg_to_ev = PCONSTANTS(7);

%vtr_distr = get_v_tr_distr(f0, 2, plist, 60);
%test = vtr_distr{10};
%plot(test(:,1), test(:,2));
%n_in_cells = get_n_in_cells(f0,plist);
%plot(n_in_cells);
t_start = plist(16);
t_end = plist(17);
i = 1;
test_pos = 110;
savedir = 'D:\Dmitry Odzerikho\heat_results\09_06_10';
savedata( plist, savedir, 'plist_current', 0);

current_pool = matlabpool('size');
if (current_pool > 0)
    matlabpool close
end
for t = t_start:dt:t_end
    pool = get_matlabpool('pool.dat');
    current_pool = matlabpool('size');
    if (pool ~= current_pool)
        if (current_pool ~= 0)
            matlabpool close
        end
        matlabpool('open', pool);
    end
    complete = 100*(t-t_start)/(t_end-t_start);
    s=sprintf('Progress: %f%%', (complete))
    f1 = simple_iter_over_distr(PCOPY, f0, t, t+dt);
    f0 = clean_distribution(f1, plist);
    vtr_distr = get_v_tr_distr(f0, 1, plist, 50);
    vtr_test = vtr_distr{test_pos};
    T_mean = get_T_mean(vtr_distr);
    savedata( T_mean, savedir, 'T_mean', i);
    savedata( vtr_distr, savedir, 'vtr_distr', i);
    savedata( f0, savedir, 'full_distr', i);
    name = get_unique_name(savedir, 'T_mean_plot', i, 'jpg');
    g=plot((z_min+0.5):1:(z_max-0.5), smooth(T_mean));
    saveas(g,name,'jpg');
    
    name = get_unique_name(savedir, 'vtr_plot', i, 'jpg');
    h=plot(vtr_test(:,1), vtr_test(:,2));
    saveas(h,name,'jpg');
    
    i = i + 1;
    
    %n_in_cells = get_n_in_cells(f1,plist);
    %plot(smooth(n_in_cells,6));
    
end
matlabpool close

answer = f0;

function f1 = simple_iter_over_distr(plist, f0, t0, t1)
nz = size(f0,1);
f1 = zeros(nz,4);
global PCONSTANTS
w_ce = PCONSTANTS(4);
%z_min = plist(5);
%z_max = plist(6);
%dz = plist(7);
t0_reduced = w_ce*t0;
t1_reduced = w_ce*t1;
z = f0(:,1);
vx = f0(:,2);
vy = f0(:,3);
vz = f0(:,4)/w_ce;
parfor i = 1:nz 
    y0 = [0 0 z(i) vx(i) vy(i)  vz(i)];
    [~,y] = evaluate_single_particle(y0, t0_reduced, t1_reduced,plist);
    f1(i,:)=[y(1) y(2) y(3) w_ce*y(4)];
end


function nz = search_cell(plist,z)
    z_min = plist(5);
    z_max = plist(6);
    
    dz = plist(7);
    if ((z > z_max) || (z < z_min))
        nz = 0;
    else
        z_curr = z_min;
        nz = 1;
        while(((z <= z_curr) || (z >= z_curr+dz)))
            nz = nz + 1;
            z_curr = z_curr + dz;
        end   
    end
    
% function distr = generate_distr_in_cell(nz, plist)
% global PCONSTANTS
% m = PCONSTANTS(2);
% Tx = plist(1);
% %Ty = plist(2);
% %C =  plist(4); %%C*v_temp = v_max
% z_min = plist(5);
% %z_max = plist(6);
% dz = plist(7);
% rho = plist(8);
% sigma = ((Tx*1.6e-12)/m)^0.5;
% N = ceil((rho^(1/3))*dz);
% distr_v = normrnd(0, sigma, N, 3);
% distr_z = dz*rand(N,1)+(z_min+(nz-1)*dz);
% distr = horzcat(distr_z, distr_v);

function [f0, meanvz] = generate_initial_distr(plist)
global PCONSTANTS
m = PCONSTANTS(2);
Tx = plist(1);
%Ty = plist(2);
%C =  plist(4); %%C*v_temp = v_max
z_min = plist(5);
z_max = plist(6);
%dz = plist(7);
rho = plist(8);
sigma = ((Tx*1.6e-12)/m)^0.5;
N = ceil((rho^(1/3))*(z_max-z_min));
distr_v = normrnd(0, sigma, N, 3);
distr_z = (z_max-z_min)*rand(N,1)+z_min;
f0 = horzcat(distr_z, distr_v);
meanvz = mean(abs(f0(:,4)));


function [t,y] = evaluate_single_particle(y0, t0, t1, plist)
[t,tempy] = ode45(@(t,y) odefun(t,y,plist), [t0 t1], y0);
y=tempy(end,3:6);

function dydt=odefun(t,y, plist)
dydt=zeros(6,1);
dydt(1) = 0;
dydt(2) = 0;
dydt(3) = y(6);
dydt(4) = -y(5) - E(t,y, plist);
dydt(5) = y(4);
dydt(6) = 0;



function E_ampl = E(t,y,plist)
z = y(3);
E0=plist(9);
E_zwidth=plist(10);
E_zpos=plist(11);
E_twidth=plist(13);
E_tpos=plist(14);
if ((abs(z-E_zpos)<2.5*E_zwidth) && (abs(t-E_tpos)<2.5*E_twidth))
    E_ampl = E0*exp(-((t-E_tpos)/E_twidth)^4)*exp(-((z-E_zpos)/E_zwidth)^2)*sin(t);
else
    E_ampl = 0;
end

function distr = get_n_in_cells(f0,plist)
    global PCONSTANTS
    z_min = plist(5);
    dz = plist(7);
    nz = PCONSTANTS(5);
    distr = zeros(nz,1);
    for i = 1:1:nz
         temp = z_min+i*dz-f0(:,1);
         a=0;
         [a,~] = find((temp > 0) & (temp < dz));
         distr(i) = size(a,1);
    end
    
function f1 = clean_distribution(f0, plist)    
    global PCONSTANTS
    m = PCONSTANTS(2);
    Tx = plist(1);
    %Ty = plist(2);
    %C =  plist(4); %%C*v_temp = v_max
    z_min = plist(5);
    z_max = plist(6);
    sigma = ((Tx*1.6e-12)/m)^0.5;
    full = 1:1:size(f0,1);
    
    [out1,~] = find(f0(:,1)>z_max);
    [out2,~] = find(f0(:,1)<z_min);
    out = size([out1; out2],1);
    distr_v = normrnd(0, sigma, out, 3);
    distr_z = (z_max-z_min)*rand(out,1)+z_min;
    df = horzcat(distr_z, distr_v);
    f1 = vertcat(f0(setdiff(full, [out1; out2]),:), df);

function vtr_distr = get_v_tr_distr(f0, DZ, plist, N_vtr)
z_min = plist(5);
z_max = plist(6);
if size(DZ,1)>0
    dz = DZ(1,1);
else
    dz = plist(7);
end
Nz = ceil((z_max-z_min)/dz);
vtr_distr = cell(Nz,1);
z = f0(:,1);
vx = f0(:,2);
vy = f0(:,3);
parfor i = 1:1:Nz
    z_min_curr = z_min+(i-1)*dz;
    z_max_curr = z_min+i*dz;
    [matched,~] = find((z>=z_min_curr)&(z<=z_max_curr));
    vtr_curr = (vx(matched).^2+vy(matched).^2).^0.5;
    min_vtr_curr = min(vtr_curr);
    max_vtr_curr = max(vtr_curr);
    dvtr = (max_vtr_curr - min_vtr_curr)/N_vtr;
    vtr_distr_temp = zeros(N_vtr,2);
    j = 1;
    for v_tr = min_vtr_curr:dvtr:max_vtr_curr
        [rows,~] = find(abs(vtr_curr(:)-v_tr)<dvtr);
        if (size(rows,1) > 0) 
            vtr_distr_temp(j,1) = mean(vtr_curr(rows));
            vtr_distr_temp(j,2) = size(rows,1);
        end
        j = j + 1;
    end
    vtr_distr{i}=vtr_distr_temp;
end

function T_mean = get_T_mean(vtr_distr)
global PCONSTANTS
m = PCONSTANTS(2);
erg_to_ev = PCONSTANTS(7);
N = size(vtr_distr,1);
T_mean = zeros(N,1);
for i = 1:1:N
    vtr = vtr_distr{i};
    norm = sum(vtr(:,2));
    if (norm > 0)
        T = 0.5*m*erg_to_ev*(vtr(:,1).^2);
        T_mean(i) = sum(vtr(:,2).*T(:))/norm;
    end
end

function name = get_unique_name(dir, name_template, num, ext)
filename = [dir, '\', name_template];
[name,~] = sprintf('%s_%i.%s', filename, num, ext);

function savedata(data, dir, filename_template,filenum)
name = get_unique_name(dir, filename_template, filenum, 'mat');
save(name, 'data');

function pool = get_matlabpool(file)
    pool = load(file);

